---
title: "Main Results (IV)"
author: "Andrei Munteanu"
date: "11/01/2022"
output:
  html_document:
    fig_width: 8 
    fig_height: 8 
  # html_document:
    # toc: true 
    # toc_float: true
---
<style type="text/css">
.main-container {
  max-width: 1800px;
  margin-left: auto;
  margin-right: auto;
}
</style>

```{r options, include=FALSE, echo=FALSE, cache=FALSE}
knitr::opts_knit$set(root.dir=getwd())
knitr::opts_chunk$set(cache.lazy=FALSE,cache=TRUE,warning=FALSE,echo=TRUE,comment='',prompt=T)
options(modelsummary_format_numeric_latex='plain')
knitr::opts_template$set(
  regular=list(include=TRUE,echo=TRUE,cache=FALSE,warning=FALSE,message=FALSE),
  regular_cache=list(include=TRUE,echo=TRUE,cache=TRUE,warning=FALSE,message=FALSE),
  invisible=list(include=FALSE,echo=FALSE,cache=FALSE,warning=FALSE,results=FALSE,message=FALSE),
  muted=list(include=TRUE,echo=TRUE,cache=FALSE,warning=FALSE,results=FALSE,message=FALSE),
  latex=list(include=TRUE,echo=TRUE,cache=FALSE))

perc.rank <- function(x) {
      y<-rank(x)/length(x)
      return(y)}

perc.rank.zero <- function(x) {
      
       y<-ifelse(x==0,0,rank(x)/length(x))
      return(y)}

data_regression2<-data_regression %>%
 mutate(dec_town=relevel(dec_town,ref="6"))
```

```{r, opts.label='muted'}
#get town/school stats including dropouts/non-matched students
data_regression2<-data_regression2 %>%
  filter(an %in% c(2008:2019)) %>%
  group_by(judet_adm,liceu_repartizat) %>%
  mutate(school_mean=mean(entrance_perc,na.rm=T)) %>%
  group_by(judet_adm,liceu_repartizat,specializare_adm) %>%
  mutate(class_mean=mean(entrance_perc,na.rm=T)) 

data_regression2$entrance_perc_town<-NA
data_regression2[!is.na(data_regression2$media_la_admitere),]<-data_regression2[!is.na(data_regression2$media_la_admitere),] %>%
      group_by(an_adm,Cod_SIRUTA2_hs_adm) %>%
      mutate(entrance_perc_town=perc.rank(media_la_admitere)) %>%
      ungroup
data_regression2<-data_regression2 %>% mutate(dec_town=cut(entrance_perc_town, breaks = c(-Inf, 0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9, Inf), 
                                                       labels = c('1','2','3','4','5','6','7','8','9','10'), right = FALSE))
    
    

#add missing values for students with missing graduation exam scores
data_regression2<-data_regression2 %>%
  group_by(Cod_SIRUTA2_hs_adm) %>%
  mutate(Wages_hs_bac=ifelse(is.na(Wages_hs_bac),mean(Wages_hs_bac,na.rm=T),Wages_hs_bac),
         drop_hs_hs_adm=ifelse(is.na(drop_hs_hs_bac),mean(drop_hs_hs_bac,na.rm=T),drop_hs_hs_bac),
         Unemployment_hs_bac=ifelse(is.na(Unemployment_hs_bac),mean(Unemployment_hs_bac,na.rm=T),Unemployment_hs_bac),
         n_students_town_yr=ifelse(is.na(n_students_town_yr),mean(n_students_town_yr,na.rm=T),n_students_town_yr),
         n_hs_town=ifelse(is.na(n_hs_town),mean(n_hs_town,na.rm=T),n_hs_town)) %>%
  group_by(judet_adm) %>%
          mutate(Unemployment_hs_bac=ifelse(is.na(Unemployment_hs_bac),mean(Unemployment_hs_bac,na.rm=T),Unemployment_hs_bac),
          drop_hs_hs_bac=ifelse(is.na(drop_hs_hs_bac),mean(drop_hs_hs_bac,na.rm=T),drop_hs_hs_bac),
          Wages_hs_bac=ifelse(is.na(Wages_hs_bac),mean(Wages_hs_bac,na.rm=T),Wages_hs_bac)) %>%
  ungroup


data_regression2<-as.data.frame(data_regression2)
data_regression2$n_hs_town_group<-data_regression2$n_hs_town
data_regression2$n_hs_town_group[data_regression2$n_hs_town>=4 & data_regression2$n_hs_town<=15]<-"4-15"
data_regression2$n_hs_town_group[data_regression2$n_hs_town>15]<-"16+"
data_regression2$n_hs_town_group<-with(data_regression2, reorder(n_hs_town_group, n_hs_town))


data_regression2<-data_regression2 %>% 
  filter(!is.na(media_la_admitere))

test<-data_regression2 %>% 
  filter(is.na(entrance_perc))

data_regression2$grad_perc_02<-NA
data_regression2[!is.na(data_regression2$media_0),]<-data_regression2[!is.na(data_regression2$media_0),] %>%
      group_by(an) %>%
      mutate(grad_perc_02=perc.rank.zero(media_0)) %>%
      ungroup







```

## {.tabset}
### Tide
```{r, opts.label='muted'}
# model_tide<-feols(grad_perc_02~entrance_perc+dec_town+n_hs_town_group+dec_town*n_students_town_yr+Unemployment_hs_adm*dec_town+Wages_hs_adm*dec_town+drop_hs_hs_adm*dec_town|
#                          as.factor(an)+specializare_adm+scoala_de_provenienta+Cod_SIRUTA2_hs_adm,
#                         cluster=~Cod_SIRUTA2_hs_adm,
#                        data=data_regression2 %>% filter(school_change==F | is.na(school_change)))
# summary(model_tide)

model_tide<-feols(grad_perc_0~entrance_perc+n_hs_town_group+dec_town*n_students_town_yr+Unemployment_hs_adm*dec_town+Wages_hs_adm*dec_town+drop_hs_hs_adm*dec_town|
                         as.factor(an)+specializare_adm+scoala_de_provenienta+Cod_SIRUTA2_hs_adm,
                        cluster=~judet_adm,
                       data=data_regression2)
summary(model_tide)
```


### Track
Run regression:

```{r, opts.label='muted'}
model_track<-feols(grad_perc_0~entrance_perc+dec_town+n_hs_town_group+dec_town*n_students_town_yr+Unemployment_hs_bac*dec_town+Wages_hs_bac*dec_town+drop_hs_hs_bac*dec_town|
        as.factor(an)+specializare_adm+scoala_de_provenienta+Cod_SIRUTA2_hs_adm|
        class_mean~dec_town*n_hs_town_group,
       cluster=~judet_adm,
      data=data_regression2)
stage1_track<-summary(model_track,stage=1)
vcov_cluster_track<-stage1_track$cov.unscaled

```

Show first and second stages

```{r, opts.label='regular'}
summary(model_track)
stage1_track
```

```{r, opts.label='muted'}
model_track_ols<-feols(grad_perc_0~entrance_perc+class_mean+dec_town+n_hs_town_group+dec_town*n_students_town_yr+Unemployment_hs_bac*dec_town+Wages_hs_bac*dec_town+drop_hs_hs_bac*dec_town|
       as.factor(an)+specializare_adm+scoala_de_provenienta+Cod_SIRUTA2_hs_adm,
                        cluster=~judet_adm,
      data=data_regression2)


```

Show first and second stages

```{r, opts.label='regular'}
summary(model_track_ols)
```
```{r, opts.label='muted'}
sizes<-unique(data_regression2[order(data_regression2$n_hs_town_group),]$n_hs_town_group)

quantile<-'dec_town'
n_quantiles<-10
deciles<-c("dec_town2","dec_town3","dec_town4","dec_town5","dec_town6","dec_town7","dec_town8","dec_town9","dec_town10")
deciles<-c("dec_town1",deciles[deciles %in% colnames(vcov_cluster_track)])


fe_stage_1_track<-data.frame()
for (i in 1:length(sizes)){
  size<-sizes[i]
  print("size")
  print(size)
  for (q in deciles) {
    print("q")
    print(q)
    q<-as.numeric(gsub("dec_town","",q))
    if (q!=1 & size!=1){
      town_name <-paste0('n_hs_town_group',size)  
      quart_name <-paste0(quantile,q)
      quart_town_name <-paste0(quantile,q,':n_hs_town_group',size)
      names<-c('entrance_perc',town_name,quart_name,quart_town_name)
      #names<-c(quart_name,quart_town_name)
      vcov_temp<-vcov_cluster_track[names,names]
      
      fe_town<-summary(stage1_track)$coefficients[town_name]
      fe_quant<-summary(stage1_track)$coefficients[quart_name]
      fe_town_quant<-summary(stage1_track)$coefficients[quart_town_name]
      
      fe<-ifelse(is.na(fe_town),0,fe_town)+ifelse(is.na(fe_quant),0,fe_quant)+ifelse(is.na(fe_town_quant),0,fe_town_quant)
      
      #fe<-ifelse(is.na(fe_quant),0,fe_quant)+ifelse(is.na(fe_town_quant),0,fe_town_quant)
    } else if (q==1 & size!=1){
      town_name <-paste0('n_hs_town_group',size)  
      names<-c('entrance_perc',town_name)
      vcov_temp<-vcov_cluster_track[names,names]
      #vcov_temp<-0
      
      fe_town<-summary(stage1_track)$coefficients[town_name]
      
      fe<-ifelse(is.na(fe_town),0,fe_town)
      #fe<-0
    } else if (q!=1 & size==1){
      quart_name <-paste0(quantile,q) 
      names<-c('entrance_perc',quart_name)
      vcov_temp<-vcov_cluster_track[names,names]
      
      fe_quant<-summary(stage1_track)$coefficients[quart_name]
      
      
      fe<-ifelse(is.na(fe_quant),0,fe_quant)
    } else if (q==1 & size==1){
      fe<-0
    }
    fe_grade<-summary(stage1_track)$coefficients['entrance_perc']*
      mean(data_regression2[data_regression2$dec_town==q & data_regression2$n_hs_town_group==size & !is.na(data_regression2$entrance_perc),]$entrance_perc)
    fe<-fe+fe_grade
    
    q<-q
    size<-size
    #se<-sqrt(sum(vcov_temp))
    #sqrt(vcov_cluster_track[town_name,town_name])+sqrt(vcov_cluster_track[quart_name,quart_name])+sqrt(vcov_cluster_track[quart_town_name,quart_town_name])+
    #            2*(vcov_cluster_track[town_name,quart_name])+2*(vcov_cluster_track[quart_name,quart_town_name])+2*(vcov_cluster_track[town_name,quart_town_name])
    if (q!=1 | size!=1){
      #fe<-data.frame(fe=fe,decile=q,n_hs_towns=size,se=se)
      se<-sqrt(t(as.matrix(c(fe_grade,rep(1,dim(vcov_temp)[1]-1))))%*%vcov_temp%*%(c(fe_grade,rep(1,dim(vcov_temp)[1]-1))))
      fe<-data.frame(fe=fe,decile=q,n_hs_towns=size,se=se) 
    } else {
      #fe<-data.frame(fe=0,decile=q,n_hs_towns=size,se=0)
      se<-fe_grade*summary(stage1_track)$se['entrance_perc']
      fe<-data.frame(fe=fe,decile=q,n_hs_towns=size,se=se) 
    }
    
    #print('dim')
    #print(dim(fe)[2])
    fe_stage_1_track<-bind_rows(fe_stage_1_track,fe)
  }
}
fe_stage_1_track$fe<-fe_stage_1_track$fe-fe_stage_1_track[fe_stage_1_track$decile==1 & fe_stage_1_track$n_hs_towns==1,]$fe
```

Track 1st Stage Graph:

```{r,  opts.label='regular'}
lb<-min(fe_stage_1_track %>% mutate(lb=fe*100-se*qnorm(0.975)*100) %>% ungroup %>% select(lb))
ub<-max(fe_stage_1_track %>% mutate(ub=fe*100+se*qnorm(0.975)*100) %>% ungroup %>% select(ub))
graph_track_1<-ggplot(data=fe_stage_1_track,aes(x=decile,y=fe*100,color=n_hs_towns))+
  geom_errorbar(aes(ymin=fe*100-se*qnorm(0.975)*100, ymax=fe*100+se*qnorm(0.975)*100),  width=0.5,size=0.8,position=position_dodge(0.5))+
  geom_point(size=1,position=position_dodge(0.5))+
   theme(strip.text.x = element_text(size=18),
        axis.text.x=element_text(angle = 0, vjust=0.5,size=14),
        axis.text.y = element_text(size=14),  
        axis.title.x = element_text(size=14),
        axis.title.y = element_text(size=14),
        legend.title=element_text(size=18),
        legend.text=element_text(size=18))+
  #ylim(-0.1,0.3)+
  xlab("Own Admission Score Decile (Town-Cohort)")+
      ylab("Peer Admission Score Percentile (National-Cohort)")+
  scale_x_continuous(breaks=0:10)+
  scale_y_continuous(breaks=seq(from = -20, to = 60, by = 5),limits=c(lb,ub))+ 
  theme(legend.position = "none")
graph_track_1
```
Differences:

```{r,  opts.label='regular'}
#Bottom decile difference in peer means (across diff town sizes)
fe_stage_1_track[fe_stage_1_track$decile==1 & fe_stage_1_track$n_hs_towns==1,]$fe-fe_stage_1_track[fe_stage_1_track$decile==1 & fe_stage_1_track$n_hs_towns=='16+',]$fe
#Top decile difference in peer means (across diff town sizes)
fe_stage_1_track[fe_stage_1_track$decile==10 & fe_stage_1_track$n_hs_towns==1,]$fe-fe_stage_1_track[fe_stage_1_track$decile==10 & fe_stage_1_track$n_hs_towns=='16+',]$fe
```

Track 2nd Stage Graph

```{r,  opts.label='muted'}
name<-"fit_class_mean"
fe_stage_2_track<-fe_stage_1_track
fe_stage_2_track$se<-fe_stage_2_track$se^2*summary(model_track)$coefficients[name]^2+
  fe_stage_2_track$se^2*summary(model_track)$coefficients[name]^2+
  fe_stage_2_track$fe^2*summary(model_track)$se[name]^2
fe_stage_2_track$se<-sqrt(fe_stage_2_track$se)
fe_stage_2_track$fe<-fe_stage_2_track$fe*summary(model_track)$coefficients[name]

lb<-min(fe_stage_2_track %>% mutate(lb=fe*100-se*qnorm(0.975)*100) %>% ungroup %>% select(lb))
ub<-max(fe_stage_2_track %>% mutate(ub=fe*100+se*qnorm(0.975)*100) %>% ungroup %>% select(ub))
graph_track_2<-ggplot(data=fe_stage_2_track,aes(x=decile,y=fe*100,color=n_hs_towns))+
  geom_errorbar(aes(ymin=fe*100-se*qnorm(0.975)*100, ymax=fe*100+se*qnorm(0.975)*100),  width=0.5,size=0.8,position=position_dodge(0.5))+
  geom_point(size=1,position=position_dodge(0.5))+
      theme(strip.text.x = element_text(size=18),
        axis.text.x=element_text(angle = 0, vjust=0.5,size=14),
        axis.text.y = element_text(size=14),  
        axis.title.x = element_text(size=14),
        axis.title.y = element_text(size=14),
        legend.title=element_text(size=18),
        legend.text=element_text(size=18))+
  #ylim(-0.1,0.3)+
  xlab("Own Admission Score Decile (Town-Cohort)")+
  ylab("Own Graduation Score Percentile (National-Cohort)")+
  scale_x_continuous(breaks=0:10)+
  scale_y_continuous(breaks=seq(from = -8, to = 16, by = 2),limits=c(lb,ub))
```


```{r,  opts.label='regular'}
graph_track_2
```

Differences:

```{r,  opts.label='regular'}

#Difference at graduation between 1st decile students
fe_stage_2_track[fe_stage_2_track$decile==1 & fe_stage_2_track$n_hs_towns==1,]$fe-
  fe_stage_2_track[fe_stage_2_track$decile==1 & fe_stage_2_track$n_hs_towns=='16+',]$fe
#Difference at graduation between 10th decile students
fe_stage_2_track[fe_stage_2_track$decile==10 & fe_stage_2_track$n_hs_towns==1,]$fe-
  fe_stage_2_track[fe_stage_2_track$decile==10 & fe_stage_2_track$n_hs_towns=='16+',]$fe

#Difference at graduation between top and bottom decile students in large towns
fe_stage_2_track[fe_stage_2_track$decile==10 & fe_stage_1_track$n_hs_towns=='16+',]$fe-
  fe_stage_2_track[fe_stage_2_track$decile==1 & fe_stage_2_track$n_hs_towns=='16+',]$fe
#Difference at graduation between top and bottom decile students in small towns
fe_stage_2_track[fe_stage_2_track$decile==10 & fe_stage_1_track$n_hs_towns==2,]$fe-
  fe_stage_2_track[fe_stage_2_track$decile==1 & fe_stage_2_track$n_hs_towns==2,]$fe
```


### School

Run regression:

```{r,  opts.label='regular'}
#rm(model_track,stage1_track)
gc()
model_school<-feols(grad_perc_0~entrance_perc+dec_town+n_hs_town_group+dec_town*n_students_town_yr+Unemployment_hs_bac*dec_town+Wages_hs_bac*dec_town+drop_hs_hs_bac*dec_town|
        as.factor(an)+specializare_adm+scoala_de_provenienta+Cod_SIRUTA2_hs_adm|
        school_mean~dec_town*n_hs_town_group,
       cluster=~judet_adm,
      data=data_regression2)
stage1_school<-summary(model_school,stage=1)
vcov_cluster_school<-stage1_school$cov.unscaled

```

Show first and second stages:

```{r,  opts.label='regular'}
summary(model_school)
stage1_school
```

```{r,  opts.label='regular'}
#rm(model_track,stage1_track)
gc()
model_school_ols<-feols(grad_perc_0~entrance_perc+school_mean+dec_town+n_hs_town_group+dec_town*n_students_town_yr+Unemployment_hs_bac*dec_town+Wages_hs_bac*dec_town+drop_hs_hs_bac*dec_town|
        as.factor(an)+specializare_adm+scoala_de_provenienta+Cod_SIRUTA2_hs_adm,
        cluster=~judet_adm,
      data=data_regression2)


```

Show first and second stages:

```{r,  opts.label='regular'}
summary(model_school_ols)
```

 
```{r,  opts.label='regular'}
sizes<-unique(data_regression2[order(data_regression2$n_hs_town_group),]$n_hs_town_group)

quantile<-'dec_town'
n_quantiles<-10
deciles<-c("dec_town2","dec_town3","dec_town4","dec_town5","dec_town6","dec_town7","dec_town8","dec_town9","dec_town10")
deciles<-c("dec_town1",deciles[deciles %in% colnames(vcov_cluster_school)])


fe_stage_1_school<-data.frame()
for (i in 1:length(sizes)){
  size<-sizes[i]
  print("size")
  print(size)
  for (q in deciles) {
    print("q")
    print(q)
    q<-as.numeric(gsub("dec_town","",q))
    if (q!=1 & size!=1){
      town_name <-paste0('n_hs_town_group',size)  
      quart_name <-paste0(quantile,q)
      quart_town_name <-paste0(quantile,q,':n_hs_town_group',size)
      names<-c('entrance_perc',town_name,quart_name,quart_town_name)
      #names<-c(quart_name,quart_town_name)
      vcov_temp<-vcov_cluster_school[names,names]
      
      fe_town<-summary(stage1_school)$coefficients[town_name]
      fe_quant<-summary(stage1_school)$coefficients[quart_name]
      fe_town_quant<-summary(stage1_school)$coefficients[quart_town_name]
      
      fe<-ifelse(is.na(fe_town),0,fe_town)+ifelse(is.na(fe_quant),0,fe_quant)+ifelse(is.na(fe_town_quant),0,fe_town_quant)
      
      #fe<-ifelse(is.na(fe_quant),0,fe_quant)+ifelse(is.na(fe_town_quant),0,fe_town_quant)
    } else if (q==1 & size!=1){
      town_name <-paste0('n_hs_town_group',size)  
      names<-c('entrance_perc',town_name)
      vcov_temp<-vcov_cluster_school[names,names]
      #vcov_temp<-0
      
      fe_town<-summary(stage1_school)$coefficients[town_name]
      
      fe<-ifelse(is.na(fe_town),0,fe_town)
      #fe<-0
    } else if (q!=1 & size==1){
      quart_name <-paste0(quantile,q) 
      names<-c('entrance_perc',quart_name)
      vcov_temp<-vcov_cluster_school[names,names]
      
      fe_quant<-summary(stage1_school)$coefficients[quart_name]
      
      
      fe<-ifelse(is.na(fe_quant),0,fe_quant)
    } else if (q==1 & size==1){
      fe<-0
    }
    fe_grade<-summary(stage1_school)$coefficients['entrance_perc']*
      mean(data_regression2[data_regression2$dec_town==q & data_regression2$n_hs_town_group==size & !is.na(data_regression2$entrance_perc),]$entrance_perc)
    fe<-fe+fe_grade
    
    q<-q
    size<-size
    #se<-sqrt(sum(vcov_temp))
    #sqrt(vcov_cluster_school[town_name,town_name])+sqrt(vcov_cluster_school[quart_name,quart_name])+sqrt(vcov_cluster_school[quart_town_name,quart_town_name])+
    #            2*(vcov_cluster_school[town_name,quart_name])+2*(vcov_cluster_school[quart_name,quart_town_name])+2*(vcov_cluster_school[town_name,quart_town_name])
    if (q!=1 | size!=1){
      #fe<-data.frame(fe=fe,decile=q,n_hs_towns=size,se=se)
      se<-sqrt(t(as.matrix(c(fe_grade,rep(1,dim(vcov_temp)[1]-1))))%*%vcov_temp%*%(c(fe_grade,rep(1,dim(vcov_temp)[1]-1))))
      fe<-data.frame(fe=fe,decile=q,n_hs_towns=size,se=se) 
    } else {
      #fe<-data.frame(fe=0,decile=q,n_hs_towns=size,se=0)
      se<-fe_grade*summary(stage1_school)$se['entrance_perc']
      fe<-data.frame(fe=fe,decile=q,n_hs_towns=size,se=se) 
    }
    
    #print('dim')
    #print(dim(fe)[2])
    fe_stage_1_school<-bind_rows(fe_stage_1_school,fe)
  }
}
fe_stage_1_school$fe<-fe_stage_1_school$fe-fe_stage_1_school[fe_stage_1_school$decile==1 & fe_stage_1_school$n_hs_towns==1,]$fe
```


School 1st Stage Graph:

```{r,  opts.label='regular'}
lb<-min(fe_stage_1_school %>% mutate(lb=fe*100-se*qnorm(0.975)*100) %>% ungroup %>% select(lb))
ub<-max(fe_stage_1_school %>% mutate(ub=fe*100+se*qnorm(0.975)*100) %>% ungroup %>% select(ub))
graph_school_1<-ggplot(data=fe_stage_1_school,aes(x=decile,y=fe*100,color=n_hs_towns))+
  geom_errorbar(aes(ymin=fe*100-se*qnorm(0.975)*100, ymax=fe*100+se*qnorm(0.975)*100),  width=0.5,size=0.8,position=position_dodge(0.5))+
  geom_point(size=1,position=position_dodge(0.5))+
      theme(strip.text.x = element_text(size=18),
        axis.text.x=element_text(angle = 0, vjust=0.5,size=14),
        axis.text.y = element_text(size=14),  
        axis.title.x = element_text(size=14),
        axis.title.y = element_text(size=14),
        legend.title=element_text(size=18),
        legend.text=element_text(size=18))+
  #ylim(-0.1,0.3)+
  xlab("Own Admission Score Decile (Town-Cohort)")+
  ylab("Peer Admission Score Percentile (National-Cohort)")+
  scale_x_continuous(breaks=0:10)+
  scale_y_continuous(breaks=seq(from = -30, to = 30, by = 5),limits=c(lb,ub))+ 
  theme(legend.position = "none")
graph_school_1
```

Differences:

```{r,  opts.label='regular'}
#Bottom decile difference in peer means (across diff town sizes)
fe_stage_1_school[fe_stage_1_school$decile==1 & fe_stage_1_school$n_hs_towns==1,]$fe-fe_stage_1_school[fe_stage_1_school$decile==1 & fe_stage_1_school$n_hs_towns=='16+',]$fe
#Top decile difference in peer means (across diff town sizes)
fe_stage_1_school[fe_stage_1_school$decile==10 & fe_stage_1_school$n_hs_towns==1,]$fe-fe_stage_1_school[fe_stage_1_school$decile==10 & fe_stage_1_school$n_hs_towns=='16+',]$fe
```


School 2nd Stage Graph


```{r,  opts.label='regular'}
name<-"fit_school_mean"
fe_stage_2_school<-fe_stage_1_school
fe_stage_2_school$se<-fe_stage_2_school$se^2*summary(model_school)$coefficients[name]^2+
  fe_stage_2_school$se^2*summary(model_school)$coefficients[name]^2+
  fe_stage_2_school$fe^2*summary(model_school)$se[name]^2
fe_stage_2_school$se<-sqrt(fe_stage_2_school$se)
fe_stage_2_school$fe<-fe_stage_2_school$fe*summary(model_school)$coefficients[name]

lb<-min(fe_stage_2_school %>% mutate(lb=fe*100-se*qnorm(0.975)*100) %>% ungroup %>% select(lb))
ub<-max(fe_stage_2_school %>% mutate(ub=fe*100+se*qnorm(0.975)*100) %>% ungroup %>% select(ub))
graph_school_2<-ggplot(data=fe_stage_2_school,aes(x=decile,y=fe*100,color=n_hs_towns))+
  geom_errorbar(aes(ymin=fe*100-se*qnorm(0.975)*100, ymax=fe*100+se*qnorm(0.975)*100),  width=0.5,size=0.8,position=position_dodge(0.5))+
  geom_point(size=1,position=position_dodge(0.5))+
   # facet_wrap(~Level,ncol = 2,scales='free_y')+
      theme(strip.text.x = element_text(size=18),
        axis.text.x=element_text(angle = 0, vjust=0.5,size=14),
        axis.text.y = element_text(size=14),  
        axis.title.x = element_text(size=14),
        axis.title.y = element_text(size=14),
        legend.title=element_text(size=18),
        legend.text=element_text(size=18))+
  #ylim(-0.1,0.3)+
  xlab("Own Admission Score Decile (Town-Cohort)")+
  ylab("Own Graduation Score Percentile (National-Cohort)")+
  scale_x_continuous(breaks=0:10)+
  scale_y_continuous(breaks=seq(from = -8, to = 16, by = 2),limits=c(lb,ub))

```

```{r,  opts.label='regular'}
graph_school_2
```
Differences:

```{r,  opts.label='regular'}

#Difference at graduation between 1st decile students
fe_stage_2_school[fe_stage_2_school$decile==1 & fe_stage_2_school$n_hs_towns==1,]$fe-
  fe_stage_2_school[fe_stage_2_school$decile==1 & fe_stage_2_school$n_hs_towns=='16+',]$fe
#Difference at graduation between 10th decile students
fe_stage_2_school[fe_stage_2_school$decile==10 & fe_stage_2_school$n_hs_towns==1,]$fe-
  fe_stage_2_school[fe_stage_2_school$decile==10 & fe_stage_2_school$n_hs_towns=='16+',]$fe

#Difference at graduation between top and bottom decile students in large towns
fe_stage_2_school[fe_stage_2_school$decile==10 & fe_stage_1_school$n_hs_towns=='16+',]$fe-
  fe_stage_2_school[fe_stage_2_school$decile==1 & fe_stage_2_school$n_hs_towns=='16+',]$fe
#Difference at graduation between top and bottom decile students in small towns
fe_stage_2_school[fe_stage_2_school$decile==10 & fe_stage_1_school$n_hs_towns==2,]$fe-
  fe_stage_2_school[fe_stage_2_school$decile==1 & fe_stage_2_school$n_hs_towns==2,]$fe
```


### Latex 2nd
```{r,  opts.label='latex'}
options(modelsummary_format_numeric_latex = "plain")
f_big<-function(x) format(x, big.mark=",", scientific=FALSE, nsmall=1,digits=1)
f <- function(x) format(round(x/1000000,1) , nsmall = 1)
print(modelsummary(list(model_school_ols,model_school,model_track_ols,model_track,model_tide),
                           statistic = "std.error",
             estimate="{estimate}{stars}",
                   stars=c('^{*}'=0.1,'^{**}'=0.05,'^{***}'=0.01),
               gof_map=list(list("raw" = "nobs", "clean" = "N", "fmt" = f),
                          list("raw" = "r.squared", "clean" = "R$^2$",fmt="%.2f")),
              metrics="R2",
              output="latex",
                escape=F))





r2(model_school_ols)[2]
r2(model_school)[2]
r2(model_track_ols)[2]
r2(model_track)[2]
r2(model_tide)[2]
```

### Latex 1st
```{r,  opts.label='regular'}

summary(model_school,stage=1)
summary(model_track,stage=1)
```

```{r,  opts.label='latex'}
f_big<-function(x) format(x, big.mark=",", scientific=FALSE, nsmall=1,digits=1)
f <- function(x) format(round(x/1000000,1) , nsmall = 1)
options(modelsummary_format_numeric_latex = "plain")
print(modelsummary(list(summary(model_school,stage=1) ,
                        summary(model_track,stage=1)),
                           statistic = "std.error",
             estimate="{estimate}{stars}",
                   stars=c('^{*}'=0.1,'^{**}'=0.05,'^{***}'=0.01),
               gof_map=list(list("raw" = "nobs", "clean" = "N", "fmt" = f),
                          list("raw" = "r.squared", "clean" = "R$^2$",fmt="%.2f")),
              metrics="R2",
              output="latex",
                escape=F))


              

r2(summary(model_school,stage=1))[2]
r2(summary(model_track,stage=1))[2]
```

